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We present a first-principles approach to the calculation of the electron-phonon interaction. This approach 
solves some theoretical difficulties in the standard derivation of the electron-phonon interaction. We do not 
make a Born-Oppenheimer approximation from the outset but transform the electronic coordinates to a frame 
attached to the nuclear framework. Subsequently coupled equations are derived which connect the nuclear 
density-density correlation function to the electron Green function, the screened interaction, and the vertex. 
This set of equations is completely equivalent to the full problem and therefore higher-order effects are 
systematically included. The derived equations are further compared to those obtained from the Frohlich 
Hamiltonian. It is shown that careless use of this Hamiltonian leads to double counting but also insight is given 
why use of this Hamiltonian has led to many useful results. Finally a simple method is presented that allows 
for the inclusion of electron-phonon coupling within a density-functional context. 
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I. INTRODUCTION 

The interaction between electrons and phonons plays a 
key role in the description of many phenomena in solid-state 
physics. A wide variety of properties of solids depend on it, 
such as the resistivity of metals, the temperature dependence 
of optical spectra, and its most dramatic consequence is 
maybe that it gives rise to the phenomenon of superconduc- 
tivity. Also in molecular systems coupling between nuclear 
and electronic motion leads to variety of physical effects, 
such as rotational magnetism. 1 A more spectacular example 
is given by molecules in strong laser fields where coupling 
between electronic and nuclear motion leads to the genera- 
tion of even harmonics 2 where only odd harmonics are al- 
lowed within the Born-Oppenheimer approximation. In small 
molecules many of these phenomena can be studied using 
accurate wave-function methods. For solids methods of com- 
parable accuracy are not yet available. However, owing to 
the development of computational techniques and resources, 
a clear increase can be observed in the number of systems 
and types of properties that can be calculated from first prin- 
ciples, i.e., without the need of adjustable or phenomenologi- 
cal parameters. These calculations are usually based on 
many-body Green-function theory or density-functional 
theory. Green-function theory 3 ' 4 has been used by many 
groups to calculate quasiparticle and excitonic properties of 
solids. Density-functional theory (DFT), 5 ~ 7 on the other 
hand, has been mostly used to calculate ground-state proper- 
ties but the number of applications that calculate excited- 
state properties using time-dependent density-functional 
theory 8 " 10 (TDDFT) is increasing steadily. 11 By calculating 
Born-Oppenheimer surfaces, DFT also allows for an accurate 
determination of phonon frequencies, which are very close to 
the experimentally observed ones. In view of these develop- 
ments it would seem a natural next step to combine DFT and 
Green-function methods to calculate phenomena that depend 
on the electron-phonon coupling in a first-principles manner. 
A first step towards this goal was taken by the group of 
Gross in which DFT and Green-function methods were used 
to calculate the critical temperature of model BCS 
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superconductors. 12 Inspired by these results we aimed to de- 
rive a general scheme that allows us to calculate properties 
that depend on the electron-phonon coupling in a first- 
principles manner, using Green-function and DFT methods. 
In doing so we found several conceptual difficulties with the 
derivations of the electron-phonon coupling within the stan- 
dard literature. We subsequently found a derivation which 
resolves these difficulties and obtained a general computa- 
tional scheme to calculate the electron-phonon coupling in 
which electrons and nuclei are treated fully quantum me- 
chanically. This scheme could then be used to test the valid- 
ity of model Hamiltonians. We furthermore found a way to 
incorporate electron-phonon interactions in TDDFT linear- 
response calculations. 

The paper is divided as follows. In the second section we 
explain the difficulties associated with the standard deriva- 
tions of the electron-phonon coupling. In the third section we 
derive the form of the Hamiltonian that forms a suitable 
starting point for our derivations. In the subsequent section 
we derive the coupled equations that form the central result 
of this work. The equations are very general and are valid for 
general molecules and solids. In the fifth section we show 
how the phonons affect the effective electron-electron inter- 
action, thereby confirming results found in more phenom- 
enological ways. In the sixth section we study the validity of 
the phenomenological Frohlich Hamiltonian and show that 
careless use of this Hamiltonian leads to overscreening of the 
phonon frequencies. We subsequently show how electron- 
phonon coupling can be incorporated in TDDFT calcula- 
tions. We finally present our results and conclusions. 

II. DIFFICULTIES IN STANDARD DERIVATIONS 
OF THE ELECTRON-PHONON COUPLING 

In this section we discuss some theoretical difficulties that 
arise in the standard derivations of the electron-phonon in- 
teraction. We start out from the complete Hamiltonian of the 
electron-nuclear system 

H= f„(R) + W nn (R) + fjr) + W ee (r) + W en (R,r), (1) 
©2004 The American Physical Society 



69 115110-1 



ROBERT van LEEUWEN 



PHYSICAL REVIEW B 69, 115110 (2004) 



where 



where the Born-Oppenheimer energy surface is defined as 
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This Hamiltonian describes the interaction of N e electrons 
with N n nuclei of mass M a and charge Z a . Here T„ and T e 
are the nuclear and electronic kinetic-energy operators. The 
operators W ee , W, m , and W en describe the electron-electron, 
nuclear-nuclear, and electron-nuclear interaction, respec- 
tively. Preceding any discussion of the electron-phonon in- 
teraction two approximations are made. First of all, the 
purely electronic problem is approximately solved for fixed 
positions of the atomic nuclei. The corresponding electronic 
Hamiltonian is 



H e = T e (r) + W en (r,R ) + W ee (r), 



(7) 



where r denotes the electronic and R the nuclear positions. 
The kinetic energy is denoted by T e and the electron-electron 
and electron-nuclear interaction are denoted by W ee and 
W e „ . The corresponding Schrodinger equation is 



H e $ e (r,R Q ) = e e (R )$ e (r,R ). 



(8) 



This problem is a complicated many-body problem in itself. 
Nevertheless, in the past decade large progress has been 
made towards a first-principles solution of this problem 
mainly using many-body Green-function methods. With help 
of the so-called GW approximation 1315 good band gaps and 
spectral functions of many solids can be obtained. With a 
subsequent solution of the Bethe-Salpeter equation also good 
description of optical spectra and excitonic effects is 
possible. 11 Second, to calculate the phonons the electronic 
energy must be known for several fixed positions of the nu- 
clei from which we can calculate the Born-Oppenheimer 
(BO) energy surface. One usually writes, for a given solution 
<$> e of the clamped nuclei problem, the wave function of the 
full problem as 



16.4 



¥(r,R) = * 4 (r,R)x(R). 



(9) 



If one now optimizes ^(R) with help of the variational prin- 
ciple one obtains for \ the equation 



[T n {R) + e B0 (R)]x(R)=E X (R), 



(10) 



e B0 (R) = W nn (R) + e,(R) + (<£ J T n \<t> e ), 



(11) 



where in the last term only integration over electron coordi- 
nates is implied. 4 It should be noted that in making the varia- 
tion with respect to x a ls° mixed first-order derivative terms 
with respect to x and ^ e appear. However, for systems with 
time-reversal invariance (as we will be discussing) the wave 
functions can be chosen to be real and these terms are then 
readily seen to vanish. 17 For the calculation of phonons it is 
now assumed that the BO surface has well-defined minima 
R and that the energy surface close to these minima is well 
described by a harmonic approximation. One introduces 
deviations from equilibrium and rewrites the Born- 
Oppenheimer equation for x m terms of these coordinates. 
This gives a set of equations for coupled harmonic oscillators 
which can be diagonalized in terms of new normal coordi- 
nates Q. In these coordinates we obtain a set of independent 
oscillators, known as phonons, which have characteristic fre- 
quencies ft,- . The BO equation then attains the form 
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(12) 



where ft, is the phonon frequency corresponding to phonon 
coordinate Q, . These frequencies are obtained from the ma- 
trix of second derivatives of the energy at the minima of the 
BO surface. With the use of density-functional theory accu- 
rate energy surfaces and phonon spectra can be obtained in 
this way. 18 

The difficulties arise when we want to go beyond the 
Born-Oppenheimer approximation and calculate the effects 
of the electrons on the phonons and vice versa. To do this we 
must split up the original Hamiltonian into an electron part, a 
phonon part, and a remainder. This remainder is exactly the 
contribution that we are interested in. The usual approach is 
to simply replace the interaction between the nuclei and their 
kinetic energy by the phonon Hamiltonian and to express the 
Coulombic electron-nuclear attraction in terms of phonon co- 
ordinates. This yields 4 



H = H ph (Q) + T e (r) + W ee (r) + W en (r,Q). 



(13) 



However, with this drastic step we obtain a Hamiltonian that 
is not equivalent to the full original Hamiltonian, since the 
electronic part of the Hamiltonian was already used to deter- 
mine the phonon frequencies. We have therefore introduced 
an ill-defined amount of double counting. Another way to see 
that this Hamiltonian is not equivalent to the one we started 
with is that it is not invariant anymore under rotations and 
translations of all particles. This is because the coordinates Q 
represent internal coordinates only, as six phonon coordi- 
nates that represent the center-of-mass motion and an overall 
rotation of the system have been eliminated (they have pho- 
non frequency zero). If we regard our system as being finite, 
this means that we neglect rovibrational couplings that can 
have important effects in molecules. The derivation of the 
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Hamiltonian of Eq. (13) is therefore obviously not satisfac- 
tory. The question is therefore how to find a derivation of the 
electron-phonon coupling that does not suffer from these 
problems. An obvious way to avoid the double counting 
would be to avoid the Born-Oppenheimer approximation al- 
together and directly split the Hamiltonian in the following 
three terms: 

H = H e (r)+H n (R)+H en (r,R), (14) 

where, by adding and subtracting terms at nuclear equilib- 
rium positions R , we have defined 

H e (r) = f e (r) + W ee {r) + W f „(r,R ) + #„„(R ), (15) 

H n (R) = f „(R) + W nn (R) - W nn (Ro), (16) 

H en (r,R) = W en (r,R) - W e „(r,R () ). (17) 

The Hamiltonians H e and H n are now purely electronic and 
purely nuclear Hamiltonians and could be used to define 
"bare" electrons and phonons. The strategy would then be to 
treat the term H en in perturbation theory. However, the 
zeroth-order Hamiltonian H n for the nuclei is (apart from the 
values of the masses and charges) identical to the Hamil- 
tonian of the homogeneous electron gas and its only collec- 
tive excitation mode is a plasmon mode. This Hamiltonian 
therefore gives a completely unrealistic phonon spectrum 
which completely lacks the acoustic modes and is a bad 
starting point for perturbation theory. One could, as is done 
in several standard texts, 19 ' 20 also consider a somewhat inter- 
mediate splitting in which ions, i.e., nuclei with rigid core 
electrons attached to them, are used as basic entities. This 
does, however, not solve the problem mentioned and intro- 
duces a further approximation and an arbitrariness in the 
definition of a core electron. Some works 21 ' 22 start out from 
the full Hamiltonian, Eq. (14), to derive expressions for 
some physical quantities but subsequently assume a periodic 
symmetry, in contradiction to the full translational and rota- 
tional invariance of Hamiltonian, Eq. (14), and use an expan- 
sion of H en in phonon coordinates which again assumes a 
split-up of the Hamiltonian as in Eq. (13) leading back to the 
double-counting problem. On the other hand, we know that 
the Born-Oppenheimer approximation gives phonon spectra 
that are in excellent agreement with experimental results. 
The Born-Oppenheimer Hamiltonian, Eq. (10), should there- 
fore be a good starting point to discuss the electron-phonon 
interaction. This brings us back to the Hamiltonian of Eq. 
(13). In spite of the mentioned difficulties in its derivation, 
this model Hamiltonian has found many fruitful applications 
in the description of a wide range of effects where electron- 
phonon interaction plays a role. One usually expands W en to 
first order in Q and treats this term, which is called the 
electron-phonon interaction, in perturbation theory. Often the 
electron-phonon interaction is parametrized and the param- 
eters are determined either from experimental data or by 
physical considerations. This procedure accounts for a wide 
class of model Hamiltonians known as the Frohlich Hamil- 
tonian. Although these Hamiltonians can be very useful it is 



not clear what their range of validity is and they are unsuit- 
able for first-principles predictions. 

We conclude that from a theoretical point of view there is 
a need for a first-principles approach to the electron-phonon 
interaction that avoids the introduction of model Hamilto- 
nians. The aim of this work is therefore to provide a theoret- 
ical approach to the electron-phonon interaction that can be 
used in first-principles calculations and that does not suffer 
from the theoretical difficulties mentioned above. An impor- 
tant step in this direction has been taken by Hedin and 
Lundqvist. 14 Rather than attempting to separate the Hamil- 
tonian in an electron and a phonon part, they use the full 
Hamiltonian to derive several equations that couple the elec- 
tron Green function, the screened interaction, the vertex, and 
the nuclear density-density correlation function. By iteration 
of these equations one obtains increasingly sophisticated ap- 
proximations of all the many-body quantities involved. This 
approach also has the important theoretical advantage of al- 
lowing for an exact definition of phonons that is independent 
of the Born-Oppenheimer approximation. The phonon spec- 
trum is then defined to be the spectrum corresponding to the 
spectral function of the exact nuclear density-density corre- 
lation function and is as such an experimental observable. 
The work of Hedin and Lundqvist has, however, two draw- 
backs. First of all, they derive their equations for classical 
nuclei, described by variables for which, at a certain point in 
the derivation, quantum-mechanical commutation relations 
must be used. Second, they start out from a Hamiltonian that 
has full translational and rotational symmetry and therefore 
the reduced quantities, such as the Green function, do not 
reflect the crystal symmetry. In this work we will remove 
these two drawbacks by dealing with quantum-mechanical 
nuclei from the outset and by referring the electronic coordi- 
nates to a body-fixed coordinate frame. 

III. THE TRANSFORMED HAMILTONIAN 

Let us start with a general remark on the full Hamiltonian 
of all electrons and nuclei of Eq. (1). This Hamiltonian is 
invariant under translations and rotations of all particles. 
This means that the ground-state wave function transforms 
under a representation of the translation and rotation group. 
Together with the inversion symmetry this implies that all 
one-body quantities such as the electron density are constant 
and that two-particle correlation functions, such as the elec- 
tron Green function, only depend on the distance between 
their arguments. This is obviously not a convenient starting 
point to describe a periodic solid. The solution of the prob- 
lem is obvious: we have to transform to a coordinate system 
that reflects the internal properties of the system. We will 
therefore carry out a coordinate transformation in which the 
electronic positions will be referred to a coordinate system 
attached to the nuclear framework. 23 In doing so we will 
assume that our system is finite but arbitrarily large so that 
our approach will be generally valid for molecules and sol- 
ids. We first define the center of mass of the nuclei: 

^CMN = Ti 2 ^A. (18) 

™ nuc a 
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where M„„= 2 "M „ is the total mass of all nuclei. The new 

it it L qf IX 

electronic coordinates are then defined to be 

r\ = K(a,p,y){r-R CMN ), (19) 

where 1Z is a rotation matrix and (a,/3,y) are Euler angles 
that specify the directions of the axes of the new coordinate 
system, which will be called the body-fixed frame (for de- 
tails see Appendix A). The Euler angles are functions of the 
nuclear coordinates. The way they depend on the nuclear 
coordinates depends on the choice of our coordinate trans- 
formation. One could for instance choose the angles in such 
a way that the nuclear inertia tensor in terms of the rotated 
nuclear variables 

R!=5Z(a I Ay)(K r 8 clfJf ) (20) 

becomes diagonal, where the nuclear inertia tensor is given 
by 

a 

Note that if we write out a particle coordinate we use the 
particle label as a superscript, i.e., R a =(R" ,R" ,R"). This 
way of determining the Euler angles is a common choice in 
nuclear physics. 24 26 However, for the description of 
phonons it is more appropriate to make a choice which mini- 
mizes the coupling between rotational and vibrational mo- 
tion. This is most conveniently done by using the so-called 
Eckart conditions 27 ' 28 that are commonly used in molecular 
physics to decouple nuclear and electronic motion. 29-35 
However, since they are rarely used in solid-state physics we 
will give a brief description of these conditions. Let R a be 
the nuclear positions that minimize the total energy within 
the Born-Oppenheimer approximation. These quantities will 
be used as parameters in our coordinate transformation. We 
choose them in such a way that 

2 M a R OM = 0, (22) 

a 

i.e., we refer these positions to their center of mass. We 
further choose these quantities in such a way that the inertia 
tensor I pq (R ) of Eq. (21) evaluated for these equilibrium 
positions becomes diagonal. Then we define the Euler angles 
as function of the nuclear coordinates R a by the following 
implicit equation: 

= 2 M„R , a XR:. (23) 

a 

It is important to realize that the numbers R a are not vari- 
ables but just conveniently chosen parameters in a coordinate 
transformation. The Coulomb potential in the electron- 
nuclear interaction now acquires the form 



1 _ 1 
|r-Rr|R CMJV +^(a,yS, r )-V-R] 

1 1 

\n(a,(3,y)(R-R CMN )-r'\ |R'-r'f 

(24) 

where one has to keep in mind that the Euler angles ( a,/3, y) 
are functions of the nuclear coordinates (R t , . . . ,R N ). The 
electron-nuclear interaction in the body-fixed frame is now 
invariant under translations and rotations of the nuclear co- 
ordinates. This is readily seen. First of all the quantities R a 
— Rcmn are invariant under a translation R^Ra + a and 
hence the Euler angles are invariant as well. Let us therefore 
consider rotations. Suppose we have a set of Euler angles 
corresponding to nuclear positions R a . Let us now apply a 
rotation O to these coordinates, i.e., we have new coordinates 
R a = OR a . For these rotated coordinates we have new Euler 
angles (5,/?, y) determined by the Eckart conditions 

N 

= 2 M^ a Kn(a,p,y)0(R a -R CMN ). (25) 

a 

Now since 1Z times O is again a rotation, it can be param- 
etrized by Euler angles, i.e., we can write 

K(a,p~,y)0 = K(a,]3,y) (26) 

for some Euler angles (a,/3,y). Now since the Eckart con- 
ditions determine the Euler angles uniquely we must have 
{a,P,y) = {a,/3,y) and we find 

K( a, p, y) (R,- - R CMN ) = K{a,p, y)(R, - R CMN ) . (27) 

We can therefore conclude that the coordinates R' t are so- 
called internal or shape coordinates that are invariant under 
rotations and translations of the nuclear framework, i.e., they 
satisfy 

R;(OR+a) = Rl(R). (28) 

(A very elegant discussion of such coordinates is given in 
Ref. 32.) Therefore the potential that the electrons in the 
body-fixed frame experience from the nuclei is invariant un- 
der rotation and translation of the nuclear coordinates. This 
is, of course, exactly the purpose of a body-fixed frame. Let 
us now turn to the other terms in Hamiltonian. One can also 
see that the electron-electron repulsion and the electronic ki- 
netic energy retain the same form in the primed as in the 
unprimed coordinates. This is simply because 1Z represents a 
rotation in which the Euler angles only depend on the 
nuclear coordinates and are independent of the electronic co- 
ordinates. However, for the same reason extra terms will 
appear in the nuclear kinetic energy. These terms have a 
physical origin. If all nuclei vibrate around their equilibrium 
position then also the axes of the body-fixed frame will vi- 
brate and therefore we are viewing the electrons from a mov- 
ing frame in which fictitious or Coriolis forces appear. In a 
diatomic molecule, for instance, this means that there is a 
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coupling between rotational and vibrational modes. For a 
solid we will show in Appendix C that such rovibrational 
terms are usually vanishingly small. The transformed nuclear 
kinetic terms become 



a 2M a ' 



where we defined 



N ' dr'- 



(29) 



(30) 



We define the mass-polarization and Coriolis terms now by 

TMPC = T n ~T n . (31) 

We have now completely specified our coordinate system 
and Hamiltonian. In order not to overcrowd our formulas 
with superscripts we will from now on drop all primes from 
the electronic variables. In the new coordinate system the 
Hamiltonian is 



where 



H = H n + H e +T MPC +W en , 



N " 1 „ 1^ 1 



2 1 2&j \r- rj \ 



II Vn 1 n V V 

,'. _ v a , 1 V a P 

n„- Z, - wvr + o 2j_ 



2M a 2 th |R a -R,6 



N„ 2V, 



i j=i \K(a,p,y)(R-R CMN )-rj\ 



(32) 



(33) 



(34) 



(35) 



The term f MPC =f MP +T c is the sum of a Coriolis part T c 
and mass-polarization part T MP which have the form (see 
Appendix B) 

3 3 

T M P=Z] (fll + fl r )P e .r+ X Pr S Pe,r P e,s, (36) 
r=\ rs=\ 



C it will be demonstrated that in that case these terms are 
inversely proportional to the diagonal elements I qci (R ) of 
the inertia tensor and therefore very small for large systems. 
This is a direct consequence of the Eckart conditions and this 
is exactly why these conditions are so suitable to define the 
transformation to the body-fixed frame. By neglect of the 
mass-polarization and Coriolis terms we have now obtained 
a Hamiltonian of the same form as the original Hamiltonian 
of Eq. (1) with the exception of the electron-nuclear attrac- 
tion term which appears in the form of Eq. (35) and breaks 
the full translational and rotational symmetry. This Hamil- 
tonian will now be used as a basis of our derivations. 

IV. DERIVATION OF SELF-CONSISTENT EQUATIONS 

We will start our derivations from the Hamiltonian de- 
rived in the preceding section. For the moment we will ne- 
glect the Coriolis and mass-polarization terms. The justifica- 
tion of this for solids is, as mentioned before, explained in 
more detail in Appendix C. In the case of molecular systems 
these extra terms are often incorporated later using perturba- 
tion theory. 33 With this approximation the transformed 
Hamiltonian in second quantization is given by 



H=T„ + T„ + W, 



+ W +W 



(38) 



where 



dx i^ 1 "(x)VV(x), 



(39) 



■-j t/xrfxV + (x)f'"(x')i/'(x')iA(x)w(r,r'), 



(40) 



W e „= dxdV ifjUx)^(x)f(R r ■ -R N )2 



|R,'-r| 

(41) 

where x=(r, cr) denotes a space-spin coordinate and R,' is a 
function of the coordinates R, as given by Eq. (20). The 
electronic field operators satisfy the usual anticommutation 
relations. 419 We further defined the nuclear density matrix 



7c=2 (v f r +v r )L eir + Z, ot rs L e>r L e ^ 



(37) 



f{R v --R N )= 2 j>\(R 1 a 1 )---ft N (R N a N J 



where P e and L e are the total electronic momentum and an- 
gular momentum. The quantities /ll, v, a, j3 are functions of 
the nuclear coordinates and are further specified in Appendix 
B where it is shown that /ll and /3 are inversely proportional 
to the total nuclear mass. The form of these equations is 
independent of the way the body-fixed frame is chosen. So 
far our derivations are valid for any system of electrons and 
nuclei, i.e., ranging from small molecules to solids. Let us 
now specify that we are dealing with solids. Obviously then 
the mass-polarization terms can be neglected since they are 
inversely proportional to the total nuclear mass. If we use the 
Eckart conditions in the specification of the Euler angles then 
also the Coriolis terms in T c will be negligible. In Appendix 



X0v(R„ o-v)---^i(RiO-i), (42) 



ll n n 



where <^(Rcr) and <^,(Rcr) are nuclear creation and annihi- 
lation operators for nucleus ; and where we summed over all 
nuclear spin variables cr, . We further defined 



dV=dR, dR 



(43) 



The reason that in the electron-nuclear interaction W en a den- 
sity matrix appears is a consequence of our transformation to 
the body-fixed frame which makes, as mentioned before, R,' 
a function of all the coordinates R, . The operators T n and 
W„„ are not written out in second-quantization form here. 



115110-5 



ROBERT van LEEUWEN 



PHYSICAL REVIEW B 69, 115110 (2004) 



The only property of these operators that we will need in this 
section is that they commute with any operator that depends 
only on electronic coordinates. Note that we have not defined 
the commutation relations for the nuclear creation and anni- 
hilation operators <j)\ and . The reason is that this does not 
only depend on the type of nucleus, i.e., either bosonic or 
fermionic, but also on the way the body-fixed frame is de- 
fined. For instance, if we choose the diagonalization of the 
inertia tensor to define the body-fixed frame then the Euler 
angles are defined by a constraint that is invariant under per- 
mutations of particles of the same type. In that case <j>\ and 
(fit will have either bosonic or fermionic commutation rela- 
tions. However, if we define the body-fixed frame using the 
Eckart conditions the commutation relations for <j>] and </>, 
will be more complicated. The reason is that the Eckart con- 
ditions are not invariant under permutations of particles of 
the same type as a consequence of the introduction of equi- 
librium positions R, . The true permutational symmetry of 
the system is then masked by the choice of our coordinate 
system, but is of course not changed. Luckily we will not 
have to use the commutation relations of the nuclear creation 
and annihilation operators in this section and our results will 
be valid for any choice of body-fixed frame (apart from the 
question whether or not the Coriolis terms are negligible). 
From a more physical point of view one might argue that in 
a solid in equilibrium the exchange probability of two nuclei 
is very small so that we may regard them as distinguishable 
particles and thereby simplify the mentioned problem. How- 
ever, since in our analysis there is no absolute need for such 
approximations we refrain from doing so. For further discus- 
sions on this point for molecules we refer to Refs. 34 and 35. 
We now define a potential operator by 



p(r) = n(r)-N(r). 



(47) 



V„(r) = dVf(Ri---Hy)2 

J " i=1 



" IR'- 



(44) 



The expectation value of this operator is the Coulomb poten- 
tial due to the nuclei felt by the electrons in the body-fixed 
frame. We further define a nuclear density operator by 



With the above definitions the electron-nuclear interaction 
can be simplified to 



W = - 

Tf on 



J r-R 



(48) 



If we evaluate the expectation value of this operator in the 
mean-field sense, we see that the electron cloud now inter- 
acts with a smeared out nuclear charge, rather than S peaks 
as in the clamped nuclei approximation. This reflects the 
quantum treatment of the nuclei. 

Now we will follow the derivation of the Hedin 
equations 13 along the lines of Hedin and Lundqvist. 14 The 
strategy is to obtain self-consistent equations for well- 
defined objects such as the electron Green function, the 
screened interaction, and the vertex function. The advantage 
of this approach is that it does not depend on any perturba- 
tion expansion of the Hamiltonian. Approximations are made 
in the final self-consistent equations that contain physical 
quantities such as dressed electron and phonon propagators. 
The self-consistent equations are derived by the functional 
differentiation method. For this purpose we define an auxil- 
iary external field <p(rf) coupling to the total charge p. Our 
Hamiltonian therefore becomes 



H=T„ + W m ,- l -\ 



+ - J Jxrfx'f , "(x)f , "(x')(A(x')(/ f (x)w(r,r') 

, , n(r)JV(R) f , 
I d 3 rd 3 R , „, + d 3 rp(r)ip(rt). (49) 



r-R 



This is the Hamiltonian we will use for all our derivations. 
First we derive the equation of motion of the electronic field 
operator tf/ in the Heisenberg picture (we use the same nota- 
tion as in the paper by Hedin 13 ), 



JV(r)=-^V 2 y„(r)= fdvf(R r --R w )E Z^r-R?). 

47T J 1 i= \ 

(45) 

The expectation value of this operator is the nuclear charge 
distribution in body-fixed frame coordinates that gives rise to 
the Hartree potential felt by the electrons. This is a smeared 
out density of nuclei, as opposed to the 8 peaks that would 
arise from a clamped nuclei approximation. One clearly sees 
from Eq. (45) that if we have the full density matrix T we 
can rotate the nuclei to the body-fixed frame after calculation 
of the density matrix, rather than doing the coordinate trans- 
formation in the Hamiltonian. The electronic density is as 
usual given by the expectation value of the density operator 



4, H (xt) = V(-T ,t)ijf(x)V(t,-T Q ), 



where — T is an initial time and 



V(f 2 ,fi) = rexp -i H{t)dt 



J i. 



(50) 



(51) 



is the evolution operator [we use the time ordering since 
H{t) contains the explicitly time-dependent external field 
<p(r?)]- The equation of motion of if/ H follows from 



id t iff H {xt) = [if; H (xt),H(t)l 
Working out the commutator yields 



(52) 



«(r) = 2 f'Xx^x). 

(T 

We further define a total density by 



(46) 



id,ip H {xt)- 



(53) 
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where we used that tp H commutes with all nuclear operators. 
If we define the time-ordered product as usual, 



nfe(»)4(x'f')]=»(f-f')W«)ft(x'/') 



-0(f-f)<AMx'f')<A„(xO, 

(54) 

where 8 is the Heaviside function, then we obtain using the 
commutation relations of the field operators and d9{t)ldt 
= 8{t) 

Ud,+ iv 2 - P (rt) W&k*)^*'*')] 



■iS(x-x')S(t-t 



') + f rfV— *— 
J r-r' 



XT[p B (r'f)fe(»)«x'0] 



(55) 



The equation of motion for ifi H immediately yields an equa- 
tion of motion for the electron Green function defined as 

G(xt,x't') 

^lUi-ToJuMTo-T^TlMxQ^ix't')]^) 
(y\U(-T ,T )V(T ,-T )\V) 

Here U is the evolution operator in absence of the external 
field <p. This means that for <p = Q the definition of the Green 
function reduces to the usual one. Note that we here work 
with time-ordered Green functions in the zero-temperature 
formalism. If one would be interested in finite temperature or 
nonequilibrium phenomena our derivations can be readily 
extended by use of the Matsubara technique 3 or Keldysh 
Green functions. 36,37 Now the functional derivative of V with 
respect to <p is 



SV(tj') 
<5(p(ri?i) 



/sgn(f-f')y(f,f 1 )p(r 1 )V(f 1 ,f') (56) 



if 1 1 is inside the time interval determined by t and t ' , oth- 
erwise SVI S(p is zero. Using this expression we can readily 
prove the following equation: 

V(T ,-T )r[fe(l)4(2)] 



8<p(3) 



= - iV(T , - T ) T[p H (3 ) M 1 ) 4(2)], (57) 

where we used the short notation ; = (x, ,f 1 ). If the spin vari- 
able is left out we write ; = (r,-,f ; ). With this equation we 
obtain the following equation for the electron Green func- 
tion: 



5(1-2) + / c/3w(l + ,3) 



SG (1,2) 
8<p(3) ' 



(58) 



where l + = (r 1 ,f 1 + A) with A a positive infinitesimal and 

w(T,2) = S(t l — t2)/\r l — r 2 \- The potential V that appears in 
the equation above is 



F(l) = <p(l)+ J d3w(l,3)(p H (3)), 



(59) 



where we used the notation 



: = (V\U(-T ,T )V(T Q ,-T )A\y) 
{ } (V\U(-T ,T Q )V(T -T )\V) ■ 

The potential V( 1) therefore corresponds to the external po- 
tential and the Hartree potential due to the electronic and 
nuclear charge distributions. Now we define the self-energy 
operator 2 by the equation 



^ tl +^v 2 -v(i) 



= S(l-2)+j d3Z(l,3)G(3,2). (61) 
From the definition of the inverse Green function 

J3G(l,3)G _1 (3,2)=5(l-2) 
follows the identity 

<5G(1,2) f 5G _1 (4,5) 



(62) 



d4d5G(l,4)- 
5<p(3) J 8<p(3) 

We therefore see that we can write 2 as 



-G(5,2). (63) 



2(1,2)=-; rf3^4w(l + ,3)G(l,4) 



8G~ l (4,2) 
S<p(3) 



(64) 



In the next step we define a dielectric function e and a 
screened interaction W by 



, _ _ SV(1) _ _ f - - - S( P (3)) 
6~ 1 (1,2)= — — = S(l-2)+ d3w(l,3)- 



Scp(2) 



8<p(2) 



W(l,2)= d3w(l,3)e-'(2,3). 



(65) 
(66) 



The dielectric function measures changes in the effective po- 
tential due to charge changes induced by the external field. 
Note that it contains both changes in the electronic and 
nuclear charge densities. We wish to study them separately. 
The electronic polarization is defined as the electronic charge 
response due to the effective field 



P e (h2) 



S(n(D) 
SV(2) 



(67) 



The electronic charge can also be calculated directly from 
the Green function 
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<«(3)>=-i2 G(3,3 + ). 



(68) 



We therefore find 



P e (l,2)=-i2 



5G(1,1 + ) 
8V(2) 



^ f «5G _i (3,4) 
■iZ d3d4G(\,3)G(4,\) — 

o-i J <5y(2) 



If we define the vertex function T by 

<SG _1 (3,4) 



we obtain 



T(34;2) = 



P e (l,2)=-i2 d3d4G(l,3)G(4,l)r(34;2). (71) 

This is the first one of Hedin's equations. A second one fol- 
lows immediately from the definition of the self-energy and 
the screened interaction, 



8V(2) 



(69) 



(70) 



2(1,2)=-/ t/3£/4£/5w(l + ,3)G(l,4) 



8G~\A,2) 8V(5) 



SV(5) 8<p(3) 
= ij rf4,i5G(l,4)W(r + ,5)r(42;5). (72) 

A third equation follows directly from the definition of the 
vertex function, 



r(12;3) = <5(l-2)<5(l-3) + 



52(1,2) 



SV(3) 



<5(l-2)<5(l-3)+ d4d5 



.52(1,2) 5G(4,5) 
'<SG(4,5) SV(3) 



8(1-2)8(1-3) 

/ 



52(1,2) 

+ | d4d5d6dl -, G(4,6)G(7,5)r(67;3). 



8G(4,5) 



(73) 



This yields an integral equation for the vertex function. In 
order to proceed we investigate the relation between W and 
P „ . From the definition of W we see that 



-- f____ 5(p(4)) _ _ 
W(l,2)=w(l,2)+ d3d4w(l,3) _ w(4,2). 

J 8(p(3) 

(74) 

If we split (p) in an electronic and a nuclear part and use the 
chain rule we obtain 



W(l,2) = w(l,2) 



_ _ _ _ 8{n(4)\ 8V(5) _ _ 
+ I d3d4d5w(l,3)— — zr~ — w(4,2) 



/ 
/ 



8V(5) 8cp(3) 

_ _ _ _ 8{N(4)) _ _ 
d3d4w(l,3)— — Z r L w(4,2) 



8<p(3) 



w(l,2)+ W(l,5)P e (4,5)w(4,2) 



- d3d4w(l,3)— — ^w(4,2). (75) 
J 8<p(3) 



This looks like the usual definition of the screened interac- 
tion except for a part induced by the nuclear charge density. 
We analyze this last term of the equation: 



S(N(4)) 
S<p(3) 



i(T[AN H (4)Ap H (3)]), 



(76) 



where AA =A — (A) defines the fluctuation of an operator A. 
If we split the total charge operator in its electronic and 
nuclear part we see that 



8{N(4)) 
8<p(3) 



i(T[AN H (4)AN H (3)]) 



+ «(T[AiV H (4)An H (3)]>. 



(77) 



The first term in the last equation represents the nuclear 
density-density correlation function 



D(l,2) = -i(T[AN H (l)AN H (2)]). 



(78) 



It is exactly this response function that interests us. It repre- 
sents a correction to the electronic screened interaction due 
to nuclear density fluctuations. We will later see that for a 
solid its spectral function represents the lattice vibrations. We 
therefore want to include this function as a variable in the 
Hedin equations. This can be done with a method introduced 
by Baym. 38 We consider an extra term 



H 2 =~\ 



dRN(R)J(R,t) 



(79) 



in the Hamiltonian that couples only to the nuclear density. 
We then see that 



5(P(D) 
(5/(2) 



■i(T[AN H (2)Ap H (l)])- 



8(N(2)) 
8<p(T) 



(80) 



Furthermore we have 
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S(p(l)) _ S(N(T)) | 
SJ{2) 8J{2) 8J{2) 

Mn(l)) 8(V{3)) 



J SV(3 



SV(3) SJ(2) 



■D{\,2) 



f - - - - - - <5<p(4)) 

+ d3d4P e (l,3)w(3,4) ^rhrr - ( 81 ) 



SJ(2) 



Solution of this integral equation yields 

8{N(2)) S{p(\)) 



8<p{\) SJ(2) 



[(l-P e w)- l D] n . (82) 



If we now insert this expression into Eq. (75) we find for W 
the integral equation 

W= w + WP e w + w(l- P e w)~ l Dw. (83) 
If we solve this equation for W we find 



W(l,2) = W e (l,2)+ J d3d4W e (l,3)D(3,4)W e (4,2), 

(84) 

where we defined the electronic part W e of the screened in- 
teraction as 



W e = w(l-P e wY 



(85) 



This completes our derivation of the Hedin equations. Let us 
summarize our results in the following set of Hedin equa- 
tions: 



2(1,2) = i </3</4G(l,3)W(l + ,4)r(32;4), (86) 



W(l,2) = W e (l,2)+ J d3d4W e (l,3)D(3,4)W e (4,2) 
W e = w{\-P e w)-\ 



(87) 
(88) 



P e (l,2)=-/X </3</4G(l,3)G(4,l)r(34;2), (89) 



r(12;3) = <5(l-2)<5(l-3) 



+ | d4d5d6dl - { G(4,6)G(7,5)r(67;3). 



dG{4,5) 



(90) 



If we put D = in these equations we obtain the usual Hedin 
equations of the rigid lattice and iteration of these equations 
leads to terms that are conveniently interpreted in terms of 
Feynman diagrams. As this is described clearly in the papers 
of Hedin and Lundqvist we will not carry out such an analy- 
sis here. The coupled equations derived here describe a gen- 
eral quantum system of electrons and nuclei and can there- 
fore be used, as we will do below, to judge the range of 
validity of approximate Hamiltonians. To illustrate the gen- 



W 



VWWWWV + A/WW\A^zT)wWWV 



W e = wwwwvv = 



+ 







SO 





FIG. 1. Diagrammatic representation of the Hedin equations. 

eral structure of these equations we have displayed them in 
Fig. 1. In this figure the interactions w, W e , and W are rep- 
resented by a wiggly, zigzag, and doubly wiggly line, respec- 
tively. The Green functions are represented with black lines 
with arrows and the bare vertex with a dot. Note that we do 
not have a determining equation for D itself. For this we 
would have to study the equation of motion of the nuclear 
creation and destruction operators which, however, leads to 
rather complicated expressions. The reason is that such an 
equation of motion will lead to taking commutators with the 
A r „-body operator F . The appearance of such Af„-body terms 
is not surprising. After all, the commonly used Born- 
Oppenheimer potential of Eq. (10) is also a Af,,-body quan- 
tity. These Af„-body terms make it difficult to find practical 
equations for the nuclear density-density correlation function 
D. However, in view of the quality of the Born-Oppenheimer 
phonons, we can expect that in practice we can get a good 
approximation for D from the Born-Oppenheimer wave 
functions. This approximation can then be inserted in the 
Hedin equations above and be iterated to obtain self- 
consistent approximations for the electronic Green function. 
The exact structure and the Born-Oppenheimer form of the 
nuclear density-density correlation function D will be studied 
in the following section. 

V. THE PHONON-INDUCED INTERACTION 
BETWEEN THE ELECTRONS 

In this section we will analyze the phonon-induced inter- 
action between the electrons in more detail. Such an analysis 
provides us with more insight into the structure of the 
nuclear-nuclear correlation function and its possible approxi- 
mations. It will also enable us to make a connection with the 
phenomenological Frohlich Hamiltonian and to judge its va- 
lidity. 

In Eq. (84) of the preceding section we saw that we can 
write the screened interaction between the electrons as the 
sum of two terms, 



W(l,2) = W e (l,2) + W / ,, 1 (l,2), 



(91) 
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where we define 



W ph (l,2)= J d3d4W e (l,3)D{3,4)W e (4,2). (92) 

The first term W e in Eq. (91) represents a purely electronic 
screening. This is the only term when we restrict ourselves to 
the approximation of clamped nuclei. The second term W ph 
describes an additional screening due to the motion of the 
nuclei and contains all the information on the electron- 
phonon interaction. The structure of this term becomes more 
transparent if we introduce the operator 



AV H (rt)=- dR 



AN H (Rt) 
|r-R| 



and define the electronic dielectric function e e by 



W e {h2)= | d3w(1.3)ee (2,3). 
We can then write 



(93) 



(94) 



s V / F S *( r i)^(r 2 ) Ff(*2)Fs(r^ 

u(r l ,r 2 ,u))= km 2j t^—- TT> : — 

+ s \ u> — Ll s + iTj u> + il s — irj J 



(101) 



In a system with time-reversal invariance we can always 
choose the eigenstates of the system real. We can therefore 
assume that the functions F s are real. Then u simplifies to 



^ 2n,F ! (r 1 )F ! (r 2 ) 
u(r l ,r 2 ;co) = lim 2 "^V^T' ( 102 ) 

Using this expression we can write an expression of similar 



structure for W 



ph ' 



W ph (r l ,r 2 ;co)= lim 2 



2n s g J (r 1 ,w)g s (r 2 ,w) 



(103) 



where we defined 



W ph (l,2)= J d3d4e; l (l,3)u(3,4)e; l (4,2), (95) 

where 

iu(r l t l ,r 2 t 2 ) = {T[AV H (r 1 t l )AV H (r 2 t 2 )]). (96) 
Inserting a complete set of eigenstates yields 

iu{r l t l ,r 2 t 2 ) = 6{t l -t 2 )J J e-Ki-WFtirMrJ 

S 

+ (1^2), (97) 



where 



F I (r) = <^ i |AV(r)|M>o>, 



(98) 



and where is the ground state with ground-state energy 
E and ^ an excited state of the system with energy E s . We 
further defined Cl s = E s — E which represents an excitation 
energy. From this analysis we see that u only depends on 
time through the combination t x — 1 2 . This is simply a con- 
sequence of the time independence of the Hamiltonian and 
applies to all two-point quantities. We can therefore Fourier 
transform with respect to the relative time coordinate and 
write 



Xe e (r 4 ,r 2 ,w). 



(99) 



Using the relation 



-1 f + 3= e~' WT 
0(T)=lim— dco— — (100) 
+ 2m J- 



a> + ir] 



g s (r\,u)= dr 2 e e (r l ,r 2 ,w)F s (r 2 ). 



(104) 



Now the phonon-induced screened interaction is in a form 
suitable for analysis. Now we study the properties of the 
functions F s (r) and demonstrate that the energies O s are of 
the order of phonon energies. The main observation is that 
the operator A V(r) only consists of creation and annihilation 
operators involving nuclear coordinates. If we treat the states 
"^o and ty s in the Born-Oppenheimer approximation it 
means that A V connects only states that differ with respect to 
excitations in the nuclear part of the wave function. In the 
Born-Oppenheimer approximation these energies correspond 
to phonon excitation energies. Let us discuss this in more 
detail. 

For F s (r) we have the expression 



dR 



(^, ; | AjV(R)|^ Q ) 
|r-R| 



^ r <%|Af(R l5 ...,R W )|^ ) 
= 2 Zj] dV 



|r-Rj| 



(105) 



we can write 



where we used Eq. (45). For s=£0 this expression is nothing 
more than the matrix element ( , P J |y„(r)|\[ , o) of the poten- 
tial V„(r) of Eq. (44). This function represents the Coulomb 
potential due to the nuclei felt by the electrons in the body- 
fixed frame. We know that for the ground state of the crystal 
this function is strongly peaked around equilibrium positions 
R ; - o ■ This is a feature of the true ground state of the system 
and is independent of the BO approximation. For the func- 
tions F s (r) we can therefore do an expansion around the 
equilibrium positions and write (we suppress the arguments 
of Af) 
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_ . . £ _ f _ v <gjAf|go) 
7=1 J l r_R ojl 

+ X ZjJ £ /y(^. 5 |Af |^o) 

x(r;-r ,o-v Roj ^ j+ 

Af " z 



(106) 



where we defined 



u,,,= J rfy(^ s |Af|^ )(R;-Roj) (107) 

and used that the integral over (^ s \ AT^o) vanishes. In a 
similar way we can in Eq. (106) also include terms from 
second order and higher in the deviations from equilibrium 
which will be necessary when we want to describe anhar- 
monic effects. Until this point we did not need to use the 
Born-Oppenheimer equation. We only needed the property 
that the nuclear density is strongly peaked in the body-fixed 
frame. However, to calculate the expectation value of the 
nuclear displacement in Eq. (107) we need approximate 
forms for the excited states M^. We will do this using the 
Born-Oppenheimer approximation. However, since Eq. (107) 
is expressed in coordinates with respect to the body-fixed 
frame [see Eq. (45)] we will have to transform the Born- 
Oppenheimer Hamiltonian of Eq. (10) to the body-fixed 
frame and expand in normal coordinates. How this is done is 
explained in detail in Appendix D. We can then use the Born- 
Oppenheimer approximation for the excited-state wave func- 
tion in terms of normal coordinates Q and write 



¥ f = *,(r,Ro)x,(Q). 



(108) 



If there are N„ nuclei then there are N=3N„ — 6 normal 
coordinates. The nuclear wave function is explicitly given in 
terms of these coordinates by 



(109) 



and where the subindex ; in Q t is a multi-index / = (q,\) for 
a normal coordinate characterized by wave vector q and po- 
larization direction A.. The indices i l , . . . ,i N label the par- 
ticular excited state x** l - G -> Xs = \h' ' an( l £n I s explic- 
itly given by 



UQ i ) = H n (Qi)e- n ^' 2 , 



(110) 



where H n is a Hermite polynomial. For the ground state we 
have jfc=0. The electronic wave function <S> e in Eq. (108) is 
the ground-state wave function of the clamped nuclei Hamil- 
tonian of Eq. (8). Excited-state electronic wave functions 
need not be considered since they have no overlap with ^Pg 
in the matrix element of Eq. (107). To keep the presentation 
as simple as possible we will from now on only consider 



monoatomic lattices with one nucleus with mass M and 
charge Z per unit cell. The normal coordinate is then defined 
by the relation 4 ' 19 



1 

R; _ R(),/ =_ /= 

\N C M 



where N c is the number of unit cells per unit volume and e qX 
is the polarization vector of the phonon. By inserting Eq. 
(Ill) into Eq. (107) we find that U; s is equal to 



*j.. = -i= 2 e q ',x'e iq ' R °-' f dV(V s \Af\V )Q q , x . 

(112) 

In Appendix D it is shown that 



JV(^ 5 |Ar|^ >e,.= dQx*(Q)G f x (Q). 



(113) 



where we defined dQ = dQ l - ■ -dQ N . The latter integral is 
readily evaluated to be 



dQx7(Q)QiX (Q)= J dQi^iQdQMQi) 

= {2VL i y m (114) 

whenever ^ J =|0---l---0) with 1 at position ; and zero 
otherwise. We therefore only need to consider singly excited 
modes. If \ s now corresponds to a state in which mode q\ is 
singly excited we obtain from Eqs. (112) and (114) 



1 



E q,V 



V2JV c Mf2 q , x 

Inserting this expression into Eq. (106) we finally obtain 



(115) 



W p/l (r 1 ,r 2 ,w) = 2 D qX (u))g qX (r 1 ,w)g* x (r 2 ,w) 



and where 



2H 



q,\ 



(116) 



(117) 



We further defined 



g q x(r,co) = (2MN c n qk r m J, J dr Y e e \r,r i;( o) 



X «q,X' V 



Ri 



iq-R 0i 



(118) 



0./I 



Expression (116) represents the effective interaction between 
the electrons that plays such an important role in the theory 
of superconductivity and has been derived before by Hedin 
and Lundqvist. However, their derivation is inconsistent 
since their final expression for the nuclear density-density 
correlation function has the periodic lattice symmetry 
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whereas their starting point is a completely rotationally and 
translationally invariant system. This inconsistency is re- 
paired in the current derivation. 

VI. COMPARISON TO THE FROHLICH HAMILTONIAN 

Using the results of the preceding section we will now 
discuss the validity of the Frohlich Hamiltonian. This Hamil- 
tonian was developed by Frohlich in the 1950s (Refs. 39- 
41) and has since then been widely used in different forms to 
describe phenomena that depend on the electron-phonon 
interaction. 20 We shall show here that the phenomenological 
Hamiltonian of Eq. (13) indeed leads to overscreening of 
phonon frequencies when higher-order terms are taken into 
account. To do this we will derive also some exact coupled 
equations for the Frohlich Hamiltonian and compare them to 
the exact equations derived in the preceding sections. To 
make the discussion as simple as possible we again consider 
a monoatomic lattice of nuclei with mass M and charge Z. 
The Frohlich Hamiltonian is then given by 

H=H ph +H e +H e _ ph , (119) 
where the phonon Hamiltonian is given by 

H ph = 2 \p\/^+\^Q^Q^ (120) 

where ju designates the polarization direction and the q sum- 
mation is restricted to the Brillouin zone denoted as BZ. The 
operators P qfL and Q q/Ji satisfy the usual canonical commu- 
tation relations as well as the relations P' q = P _ qfl and 

Qq/jL = Q-qp.- Here we used a notation close to that of Refs. 
4 and 19. The electron-phonon interaction is given by 

2 I d 3 ry q/ ,(r)n(r)Q qfl , (121) 

where the q summation now extends over all wave vectors. 
In this expression the function y qAt (r) is explicitly given by 

^v (r) = z /W e " >r ^-" ( q ) q- e w > (122) 

Vq \jMN c 

where V is the volume of the unit cell and v e -„(q) is the 
Fourier transform of the electron-nuclear interaction explic- 
itly given by 

r e ~' qr AttZ 
u e _„(q)=-Zj d 3 r-^-=- — . (123) 

The function y q/l further satisfies the useful relation y qfl 
= 7-qfj.- Note that we will follow Ref. 42 by denoting every 
Hamiltonian of the form in Eq. (119) as Frohlich Hamil- 
tonian, although often this term is reserved for special 
cases 20 of this expression. This nomenclature turned out to 
be most suitable for our general discussion. We now define a 
nuclear charge density operator by 



Mr)=-^S y qfl (r)Q qfl j. (124) 
In other words 

^ « f , N(R) 

2 r q/t (r)G qM =-J d (125) 

We then define also p(r) = n(r)—N(r) as the operator for 
the total negative charge and we will introduce an external 
field tp(rt) coupling to this charge. In this way the full 
Hamiltonian is given by 

H = H ph -\j rfx^(x)V 2 <A(x) 

+ Jxdx'f , "(x)iA + (x')i/f(x')iA(x)w(r,r') 

f , , «(r)iV(R) f , 
-J d 3 rd i R-^-^-+ J d\p{r)q>{rt). (126) 

This Hamiltonian is, apart from the replacement of T„ 
+ W nn by H ph and a different interpretation of the nuclear 
density, of identical form as the Hamiltonian of Eq. (49). The 
derivation of the Hedin equations using the functional de- 
rivative method therefore proceeds completely analogously 
as in Sec. IV and we obtain the identical Hedin equations 
(86) -(90), albeit with a different interpretation of the nuclear 
density-density correlation function D which is now evalu- 
ated with the density operator of Eq. (124) as 

D(T,2) = -i(T[AN H (T)AN H {2)]). (127) 

We further define the equivalent of operator A V H (Rf) of Eq. 

(93) as 

AV fl (rf) = -J |r _ R| =2 y q/ ,(r)A<2 qM/ (f), 

(128) 

where we also defined the fluctuation operator AQ q/J . H (t) 

— (Qqa,H(t)) ■ In analogy with the preceding section we now 
consider the time-ordered expectation value of this operator 
given by 

iu(r l t l ,r 2 t 2 ) = {T[AV H (r 1 t l )AV H (r 2 t 2 )]) 
= 2 Tq^ri) 

x(r[A<2 q/t (f 1 )A<2 q> ,(/ 2 )]> r * > ,(r 2 ). 

(129) 

We see that u and D can be calculated if we have an expres- 
sion for the quantity 

idqq^ fl ,(t 1 J 2 ) = {T[AQ qfl (t 1 )AQ q , fi ,(t 2 )]) (130) 
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which we will refer to as the phonon propagator. We can 
derive an equation of motion for this quantity if we consider 
the equation of motion of Q qfjL (t). In order to do this we will 
also add an external field to the Hamiltonian of the form 



H ext (t)= X vwe^+^we^, (i3i) 



D 



/iqe5Z 



where J^ /l =J_ vl and where J q+Gfl = J qft with G a 
reciprocal-lattice vector. We then have 

d t Q^(t)=-i[Q q 4t)M] = P{ fl (t). (132) 

The second derivative is then readily evaluated from the 
commutator with P q/Ji (t) to give 



d 2 ,Q^{t ) = d,P'^{t) = - Cl^Q^t) - J rf 3 rr*,(r 

-[J-^(t) + J^m (133) 

If we take the functional derivative with respect to J q t 
(where we treat J and J* as independent) we obtain 

(^+n 2 ^)d m ,^,(t,t')=-s(t-t')s m ,d^, 



S{n H {rt)) 



(134) 



where S qq , = 1 if q and q' differ by a reciprocal-lattice vec- 
tor and zero otherwise. Since the density is given by the 
diagonal of the Green function we can write the latter equa- 
tion as 



S(t-t')S qq ,S^,-i^ I d 3 rd3d4y^(r) 



XG(rta,3)G(4,rta) 

This equation can be rewritten as 

(^ + 0, 2 q/ Jd m , tl/1 ,(t,t') 



SG' 1 (3,4) 



(135) 



.+*2 2 f 

qja cr J 



d 3 rdt 1 d3d4y* fl (r 



XG(r/cr,3)G(4,rt<r)r(34;q 1 a* 1 )£/ qiq ,, B/ ,,(* 1 ,f'), 

(136) 

where we defined the vertex function as 



r(34;qaf) ; 



<5G _1 (3,4) 
S{Q qa {t)) ' 



(137) 



•o 



+ 



-O 



+ 



<5£ 
5C7 



FIG. 2. Diagrammatic representation of the self-consistent equa- 
tions that determine the phonon propagator for the Frohlich Hamil- 
tonian. 

This is a Dyson-like equation for the phonon propagator with 
the "bubble" GGT . Now the electronic Green function sat- 
isfies 

G _1 (l»2) = Go 1 (l,2)-2(l,2)-^(l-2) 

X ( d3w(l,3)(p(3)) 
= G 1 (U)-2(l,2)-<5(l-2) 



(138) 



X d3w(l,3)<«(3))-5(l-2) 



xE r qM (r 1 )(e q<11 (fi)) 

qfi 



Therefore 



r(l2-q^t) = S(l-2)S(t 1 -t) 7qfl (r 1 ) 

,<5£(1,2) 8G(4,5) 



+ d4d5- 



S °(^ S(Q qfl (t)) 
S(l-2)S(t 1 -t)y q/jL (r 1 ) 
52(1,2) 



+ d4d5d6dl 



xT(61;qfit), 



SG(4,5) 



G(4,6)G(7,5) 



(139) 



where it should be noted that we will also regard the phonon 
propagator d to be a functional of G. 

Let us now discuss the results that we obtained. The cen- 
tral equations of this section are summarized diagrammati- 
cally in Fig. 2. The top line in this figure expresses the 
nuclear density-density correlation function in terms of the 
phonon propagator and is explicitly given by 



D(l,2)= 2 «q,a(ri)^qqV ) "'( f l' f 2) Q; qV'( r 2). 



qq fj.fi 



(140) 



where we used Eqs. (127), (124), and (130) and we defined 

V 2 



«q i c( r )=-^rq /i (r). 



(141) 
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In the figure the quantity a qjlt (r) is denoted by a small circle 
and the phonon propagator by a dotted line. In the second 
line of Fig. 2 we represent Eq. (136) in Dyson form where 
the dashed line represents the bare propagator d° that satis- 
fies Eq. (144) and is discussed in more detail below. The last 
line of the figure represents Eq. (139) in which the bare 
vertex y qM is represented by a black dot. The equations rep- 
resented in Figs. 1 and 2 completely determine all the prop- 
erties of the Frohlich Hamiltonian. Before discussing itera- 
tive solutions of these equations we first study the structure 
of the phonon propagator. The phonon propagator can be 
written in its Lehmann form as 

S 

X<¥,|A£J v |¥ > + (1~2). (142) 

If we take the phonon frequencies from a calculation of the 
Born-Oppenheimer surface and denote the corresponding 
propagator as d we find using Eq. (114) that 

+ (l<->2). (143) 
This propagator satisfies the equation of motion 



(144) 



If we Fourier transform with respect to t x — t 2 we obtain 
simply 



The function u is then given in frequency space by 



(145) 



i«(r 1 ,r 2 ,ta)= 2 T^Oi) 



qq p/i 



2n qM iy(r 1 )F q ;(r 2 ) 

/.qeBZ C0 2 -(fl q/Ji -i7j) 2 



(146) 



where 



F qX (r) = (2n q J- 1/2 2 y q+G , x (r). (147) 

G 

If we work out the latter term we have 
f>(r) = (2fl qX MAg- 1/2 

X^S 6 qX - f rf3 r ^i(q+G)-(r-r') V ^ 

Vq g J |r | 

= (2il qX M^.)- 1/2 e q x-S J d 3 r'e"»-e-'') 



X^r-r'-RoJV- 



: (2n qX M^V c )- 1/2 S e iq Roi e q x V T 



Rn 



(148) 



where we used that e q = e q+G . We see that with this approxi- 
mation we obtain for the function g qfl an identical expression 
as in Eq. (118) of the preceding section. 

Results will differ, however, as we sum the phonon propa- 
gator of the Frohlich Hamiltonian to higher order. For in- 
stance, if we start iterating the coupled equations by taking 

for r the first term on the right-hand side of Eq. (139) we 
obtain for the phonon propagator the equation 



- S(t-t')S (l(l ,S fl/jL 



qi« tr 



d 3 rd3y^[r) 



XG(rtcr,3)G(3,rta) 7 ^ a (r 3 )d qiq ^ afir (t 3 j'). 

(149) 

This equation gives d as a functional of G which can, using 
Eq. (140), subsequently be inserted into the Hedin equations 
of the Frohlich Hamiltonian to obtain higher-order approxi- 
mations for the self-energy X and thus a new F from Eq. 
(139). The first iteration in Eq. (149) amounts to summing all 
diagrams for the phonon propagator that contain an elec- 
tronic polarization bubble P e = —iGG which represents the 
electronic screening. By this dressing of the phonon propa- 
gator one obtains new screened phonon frequencies from its 
poles. However, since we know that the Born-Oppenheimer 
phonon frequencies are very close to realistic phonon fre- 
quencies, this leads to overscreening. This illustrates the 
double-counting problem discussed in the Introduction. We 
therefore obtain an important result. If we use Born- 
Oppenheimer phonon frequencies, then the exact coupled 
Hedin equations of Sec. IV and the Frohlich Hamiltonian 
yield identical results provided that we do not dress the pho- 
non line in the Frohlich model. What does this imply for the 
results obtained with the Frohlich Hamiltonian? First of all, 
many results obtained with this model were aimed at quali- 
tative rather than quantitative agreement with experiment 
and often contain adjustable parameters that can be fitted to 
experimental data. On the other hand the calculations that 
used accurate frequencies and that aimed for quantitative 
agreement with experiment were often carried out in first- 
order perturbation theory in which double counting cannot 
occur. Nevertheless, it is important to be aware of the 
double-counting problem when one studies simple cases of 
the Frohlich model that can be solved exactly. 43 The problem 
of double counting is, as explained, related to the fact that 
the phonon frequencies were determined in a way that al- 
ready used the electronic part of the Hamiltonian. Indeed 
several works introduced bare phonon frequencies on the ba- 
sis of a Hamiltonian as in Eq. (16) and one often considers 
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effective interactions between "ions" which are loosely de- 
fined as nuclei to which core electrons are rigidly attached. 
The disadvantage of the latter approach is that it is difficult 
to give it a rigorous basis since it involves the somewhat 
arbitrary definition of a core electron. Moreover the bare 
phonon frequencies one starts out with are very unphysical. 
One has to do an infinite summation to go from unrealistic 
ionic to more realistic dressed phonon frequencies. Apart 
from being a cumbersome procedure it is very unclear if the 
exact phonon frequencies (defined by the spectral function of 
the exact D) can even be obtained by such a procedure. 
Another disadvantage of the Frohlich model as compared to 
the coupled equations of Sec. IV is that anharmonic effects in 
the electron-phonon coupling are absent. By using the 
coupled Hedin equations such effects can readily be included 
by expanding Eq. (106) to higher orders in Rj— R j. The 
general point we like to make is that in first-principles cal- 
culations the direct use of the coupled Hedin equations to- 
gether with the BO approximation for the nuclear density- 
density correlation function presents a simple and general 
way of generating self-consistent approximations for elec- 
trons and phonons that are devoid from double-counting 
problems. 

VII. PHONONS IN TIME-DEPENDENT 
DENSITY-FUNCTIONAL THEORY 

Density-functional theory, as usually applied to electronic 
systems, expresses all observables as functionals of the elec- 
tron density. In the case that we deal with the coupling be- 
tween electronic and nuclear motion one needs an extension 
of electronic density-functional theory. Such an extension 
has been provided by the multicomponent density-functional 
approach of Kreibich and Gross 23,44 and has been applied to 
diatomic molecules in strong laser fields. This theory has, 
however, not yet been investigated for the case of phonons. 
On the basis of the Hedin equations of Sec. IV one can, 
however, quite easily construct phonon corrections to results 
found with electronic density-functional theory. For instance, 
the electronic density-response function is readily calculated 
in time-dependent density-functional theory from 9 ' 10 ' 45 



x(l,2)=Xo(l,2)+ J d3d4 Xo (l,3) 

X[w(3,4)+/ A£ ,(3,4)]^(4,2), (150) 

where xo is the noninteracting density-response function of 
the Kohn-Sham system and f xc is the so-called exchange- 
correlation kernel. From this function one can then readily 
construct W e = wxw and e e . Therefore by calculation of the 
phonon frequencies using the Born-Oppenheimer approxi- 
mation, which can also be done within DFT, 1 
struct the effective interaction 



we can con- 



W(l,2) = W e (l,2) + W ph (l,2) (151) 

from which the full dielectric function e(l,2) can be ob- 
tained. This allows us then to calculate the phonon broaden- 
ing of absorption spectra. The approach here is quite differ- 
ent from more ad hoc ways to calculate the electron-phonon 



interaction within a DFT framework. Other works have done 
this by calculating the change in the effective Kohn-Sham 
potential v s due to a shift of the nuclei. 46-49 This amounts to 
using v s rather than e^ 1 in the expression for W ph (l,2). No 
justification for this procedure is given and from our results 
it appears that this procedure is indeed not justified and can 
only be regarded as an unfounded but simple approximation 
to e e 1 . The results, however, show that this approximation 
works quite well. 

VIII. CONCLUSIONS 

In this paper we have derived coupled equations for 
many-body Green functions, effective interactions, and ver- 
tices for the full system of electrons and nuclei. This ap- 
proach was inspired by existing theoretical difficulties in 
standard derivations of the electron-phonon interaction. They 
involved (a) the breaking of the full rotational and transla- 
tional symmetry of the original Hamiltonian and its corre- 
sponding neglect of rovibrational couplings and (b) the 
double-counting problem due to the use of Born- 
Oppenheimer phonon frequencies. The first problem was 
solved using a coordinate transformation that refers the elec- 
tronic motion to a frame fixed to the nuclear framework. The 
second problem was solved by deriving coupled equations 
for observable quantities such as Green functions, rather than 
to try to define bare electrons and phonons and to expand 
parts of the Hamiltonian. We further showed a way to calcu- 
late the electron-phonon coupling within a density-functional 
framework and pointed out some problems in earlier work 
that calculates the electron-phonon coupling within a DFT 
context. 

We hope that the work presented here will provide a use- 
ful basis for future first-principles approaches to the calcula- 
tion of the electron-phonon interaction. Work on applications 
within a time-dependent density-functional context is in 
progress. 

APPENDIX A: DEFINITIONS AND USEFUL RELATIONS 

The Euler angles are specified by the fact that any rotation 
can be written as 



(Al) 



where 1Z Z and TZ y are rotations about the z and y axes. More 
details can be found in Refs. 29, 33, and 50. Let the column 
vectors of this rotation matrix be denoted e l5 e 2 , e 3 , i.e., 
(ej)i = Hij. These vectors satisfy 



- = e,x Wa , 



= e ; X w p , 



— = e,Xo) 

dy ' 



(A2) 



(A3) 



(A4) 



where we defined the following angular velocity vectors: 
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( 


— sin a cos y\ 




(o\ 


(Op = 




sin a sin y 







\ 


cos a 1 







(A5) 

Now all the angles depend on nuclear coordinates 
= (R k ,R k R k ) where k labels the particle. Let us further in- 
troduce the short notation d k =dldR k . Then from the equa- 
tions above we find 

^(e,), = 2 e Jqr M q a k ri , (A6) 

qr 

where is the antisymmetric Levi-Civita tensor and where 
we defined 

n*- = ft> a<r <?fa+ (0^/3+ m y>r d\y (A7) 

and where co ir is the rth component of vector w, . Since the 
Euler angles are invariant under translation of all nuclear 
coordinates we can easily derive the useful condition 



If we now define 



X n*=o. 

k=\ 



= 2 T^jfi 
jk 



then we find the relations 



e- r'fr 

t jqr'q* L ri ' 



(A8) 

(A9) 
(A10) 

(All) 
(A12) 



These relations will be useful in the following sections. 

APPENDIX B: GENERAL FORM FOR THE CORIOLIS 
AND MASS-POLARIZATION TERMS 

Here we will derive the general form of the mass- 
polarization and Coriolis terms presented in Sec. III. In that 
section we rotated the electronic coordinates to a frame that 
is attached to the nuclei but left the nuclear coordinates in- 
tact, i.e., 



r;=ft(a,/3, r )(r-R CMW ), 

r;=r, 



(Bl) 
(B2) 



The wave function "V in the old coordinates is then related to 
wave function $ in the new coordinates by 



*(rI-.T^,R{---R^) = *(r 1 ..Tv,R r 



Ra 



(B3) 



If we use 

we obtain 
#9 r?<J> 



^0 = 2 e jqr r' q Q, k Ti -Kj t 



Mi 



(B4) 



d<& dr'J 



dRl dRf i=i « = i dr'f dR k ' dR\ 
M k ^ 

1V± nuc I 



(B5) 



where we used relation (B4) and where we defined the total 
electronic momentum and angular momentum operator 



V< , d 



7=1 



N., 



7 = 1 dt'j 

The kinetic -energy operator therefore becomes 



(B6) 
(B7) 



N„ 



T' 



1 



3 



3 



2 k-'2 «^..r 



M nuc I 



(B8) 



Let us work this term out. Because of the condition of Eq. 
(A8) there are no couplings between the L e and P e operators. 
We can then split the kinetic energy as 

T' n = t n + f MP +f c , (B9) 
where the mass-polarization terms are given by 

3 3 
*W=2 (f*l + fl r )Pe,r+ 2 P ' rs P ' t ,f 'e,s , (BIO) 



where 



N„ 3 



fl 



Prs- 



P' 



2M„ 



1 N " 3 8 

2 ^aS k,^/ 



2M Z k=\ 

nuc 



2M 

nuc 



(Bll) 
(B12) 



where P ( ' is the nuclear momentum in the body-fixed frame. 
The Coriolis terms are given by 

3 3 

7c=2 (vl+v r )L tir + 2 0L rs L er L L%s , (B13) 

r— 1 rs = 1 

where we defined 



^=1 zm^ /=i 



(B14) 
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k=l £M k ,= i 

Our next task will be to provide more explicit expressions 
for v r and a rs . 

APPENDIX C: EXPLICIT FORM OF THE CORIOLIS 
TERM FOR SPECIFIC CHOICE OF THE 
BODY-FIXED FRAME 

In this section we will demonstrate that for a large system 
the Coriolis terms become vanishingly small if we specify 
the Euler angles using the Eckart conditions given in Eq. 
(23). Since 

= 2 M i R 0J xn(a,(3,y)R CMN (CI) 
the Euler angles are equivalently defined by the conditions 

= 2 M i R 0i XTl(a,p,y)R i . (C2) 

i = i 

This simply means that the angles are invariant under trans- 
lations of the original coordinates. Classically the Eckart 
conditions enforce that in the rotated body-fixed frame the 
internal nuclear angular momentum with respect to the equi- 
librium configuration R, is zero. This follows in a classical 
system directly from differentiation of the Eckart conditions 
with respect to time. 29 We now calculate the angular velocity 
vectors . An equation for these vectors can be obtained 
from a differentiation of the Eckart conditions with respect to 
R k . If we denote R' k = lZR k we obtain 

= 4 2 M neijs R n 0J Rf\ 

=2 M n e ijs R"Jj, £ sgr R';"n k rl +n sl s k „) 

njs \ qr J 

= 2 M k e ijs R k 0J TZ sl +^ M n (S iq S jr -S ir S jq )R n 0J R n q '[l k rl 

js njqr 

= 2 M k e ijs R k 0J K sl +J, M n (R n OJ Rf^jrR"ojRf^ii) 

= 4-2 /y(R")0*/, (C3) 
where the matrix J(R) is defined as 

J pq (R) = 2 M k [(R or R k )S pq -R k 0j ,R k ] (C4) 

k — 1 

and where 

4=2 M^XyT^^Ro^Xe,),.. (C5) 



Note that the matrix J \, is symmetric because of the Eckart 
conditions. Furthermore, the matrix J(Rq) is diagonal be- 
cause of the conditions posed on the R , and corresponds to 
the inertia tensor of the points R Q k . The angular momentum 
vectors Cl k , can now be obtained from 

3 

ftj, = 2 (J~%al- (C6) 
p 

Using the explicit form of a lp and the fact that the R 0j £ are 
defined in a frame with respect to the nuclear center of mass, 
we can check that condition (A8) is indeed satisfied. Let us 
now calculate the matrix a,j of Appendix B. For this we 
need to evaluate 

3 

2 «£,<-=2 (J-XsiJ-%,1, 44- (C7) 

i st i 

The last term in this equation is readily calculated to be 

2 4*4= 2 ^suv^qrMkM^Jl 1 ^ KiKi 
i uvqr i 

= 2 € suv € !qv M k M l R 0M R 0,q 

ii vq 

= S s { 2 M k M,R k lq R l Q ^ -M k M,R k 0j R' Q , s , 

(C8) 

where we used the orthogonality of the matrix 1Z. With this 
equation we obtain 

N n j 3 

a ^ = ? 2M^? ^''^ 

=2 (-/ _1 ) K . s (^ _1 ) ? r2 4*4 

si k ^M k i 

= 7 2 (r^/JRo)^ 1 )^, (C9) 

where I SS (R ) are the diagonal elements of the inertia tensor 
at the equilibrium positions R . Since this tensor is diagonal 
we can also write a uq as a matrix product (where we keep in 
mind that J is a symmetric matrix and hence / too): 

a pq =k[J- l {R")I{R Q )J-\R n )] pq . (C10) 

Since the matrix J is invariant under translations, i.e., J(R 
+ a) = 7(R), we can rewrite this as 

a pq =KJ~ l (R')nW\^')] Pq , (Cll) 

where we translated over the vector — 7ZR CMN . It remains to 
calculate a more explicit form for the operator v r . For this 
term we find 
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1 



*=1 k i= \ k =i ZM k ip 1 



1 



^ k ipst p 



P • 

(C12) 



In order that the transformation from the old to the new 
coordinates be invertible the coefficients b\j must be chosen 
in such a way that the new coordinates Q, are orthogonal to 
the conditions 



= 2 M k Rl, 



(D4) 



where we defined the vibrational angular momentum of the 
nuclei by 



o=2 m a ,r (U xr;, 



(D5) 



•"V 



(C13) 



We see that if R'=R then a pq is diagonal and inversely 
proportional to the diagonal elements of the inertia tensor 
and also v r has the same proportionality. This is therefore a 
small quantity when the system is large. It is exactly this 
property that makes the Eckart conditions useful for the 
study of small molecular vibrations. 

APPENDIX D: NORMAL COORDINATES 
IN THE BODY-FIXED FRAME 

In this section we derive the body-fixed frame expression 
for the Born-Oppenheimer Hamiltonian that is used in Sec. V 
to construct the phonon propagator. Suppose we consider a 
nuclear Hamiltonian of the form 



ff=-E KTT^l +V(R 1 ...R JV ). (Dl) 

A=l lM k k 



In our case the potential V will describe the BO surface. We 
now consider the new coordinates 



Rl = K(a,/3,y)(R-R CMN ), 



(D2) 



where the Euler angles are determined by the Eckart condi- 
tions. 

These R,' are shape coordinates in the sense that they are 
invariant under rotations and translations of the R, . They 
therefore span a (3N— 6) -dimensional space and are there- 
fore dependent. If we remove R^ and RJ^_ l we are left with 
at most 3N — 6 independent shape coordinates. We therefore 
have to select six coordinates spanning the space "orthogo- 
nal" to the space of shape coordinates. One obvious candi- 
date is the nuclear center of mass Rcmn as defined in Eq. 
(18), which is not in the shape space since it is not transla- 
tionally invariant. For the remaining coordinates we choose 
the Euler angles a,/3,y. These coordinates are translation- 
ally invariant but not rotationally invariant. We are therefore 
left with a choice of 3N — 6 shape coordinates. For this we 
take the general form 



e,=2 2 bijRf. 



(D3) 



which are the center-of-mass fixing and Eckart conditions. 
For details on this procedure we refer to Refs. 34 and 51. We 
now have a new wave function <I> related to the old one 
by 



0>(e 1 ,...,(23^_ 6 ,R CMJV ,a,yS,r) = ^(Ri, ...,»*). 



(D6) 



Then partial differentiation yields 



dW d<i> , d<b , d<i> , 
dR k da 1 d/3 ,H dy ' r 

V d ® k V d ® k 

1=1 oR CMN, I j=l d Qj 

(D7) 

where the derivatives of the Euler angles can be expressed in 
terms of the coefficients ft by inversion of Eq. (A7). This 
yields 



(?. a = sin yfl*, + cos yd^,-, 



kn cos y^k , Sln 7 n< 



sin a since 

nk t _.r»* .■„,„; r>k <r>k 



(D8) 
(D9) 



d- y=cotacos yfl^-cotasin y^l k li + D,\ i . (D10) 
We then have 

<?<!> . d<$> , d<& . J, . 

^ a+ ip^ + ^y=-%^ (D11) 

where we defined the operators 

d cos y d 



iJ l = sin y 



+ cotacos y— , (D 1 2) 
da sin a dp dy 



d sin y d d 

;7 2 = cosy— + — - cotasin y— , (D13) 

da sm a d/3 dy 



d 

dy' 



(D14) 



which represent the rovibronic angular momentum operators 



relative to the body-fixed axes. 33 Using 
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Mi 



a?*"'=2 e^'a^+nJs^-j^) (Di5) 



we further obtain 



3 3 
7=1 



The term T CMN is simply the kinetic energy of the center of 
mass. The terms T R and T vih describe, loosely speaking, the 
rotational and vibrational part of the kinetic energy. If finally 
we choose the matrix fey- in such a way that the vibrational 
kinetic energy becomes diagonal and that the potential V has 
a harmonic expansion in Q i and neglect higher-order terms 
we obtain for H vlh = T vib + V the form 



where we defined 



2 HrfPr*. (D16) 



3 AT 



M 



r = 1 m = 1 



3N-6 



„m ; X i m 



L N =-i^Z R'"'xp'". 



Collecting our results together we obtain 



(D17) 



(D18) 



3N-6 

)■ = 1 



1 d 2 

2 dQ\ 



1 



(D23) 



where the fl, are identified with the vibrational phonon fre- 
quencies. The kinetic energy of the total nuclear center of 
mass can now be separated from the Hamiltonian. If we ne- 
glect Coriolis and mass-polarization terms and evaluate fly 
at the equilibrium positions Rq then T R can be replaced by 
the rigid-rotor-type Hamiltonian 



Tr 2 a rsJ i J s 



(D24) 



3 



■»2 n k ri (J r +L N j<i>-iJ, n riP \<s> 



3 N 



M k d<& 

+ h? i2 2 w$ 

M nuc \ dK CMN,i r=\ m=\ j 



(D19) 



The transformed kinetic energy is therefore given by 



N . 3 

r=-y — — y 



+ 



M t 



M nuc \dR CMN j r =i m =i 



3 3 

*2 KiPl-'lZ n k ri (Jr + L N ,r) 
r=\ r=l 

3 AT \ 12 

il X . (D20) 



where a rs is given by expression in Eq. (B15). Using the 
Eckart conditions we then find that a rs is inversely propor- 
tional to the diagonal elements of the inertial tensor. We can 
now write the Hamiltonian as 



H— T CMN + T R + H V 



!vib . (D25) 
The eigenfunctions of this Hamiltonian are of the form 



V q <JLctiN)<Pi(a,P,y)X,(G), (D26) 

where rj q is a plane wave corresponding to the center-of- 
mass motion 



1 



'7 q (RcMiv)= ~l=e 



iq R c 



V 



(D27) 



Let us work out this expression. Because of condition, Eq. 
(A8), there are no mixed expressions between the terms in- 
volving the angular momentum J + L N and the center-of- 
mass/polarization terms with prefactor M k IM nuc . We obtain 



T ~ Tcmn^ Tr + T vib , 



(D21) 



where 



which is normalized in a volume V that will cancel out of our 
final equations. The function <p/ is an eigenfunction of the 
rigid-rotor Hamiltonian of Eq. (D24) and \ s is an eigenfunc- 
tion of the vibrational Hamiltonian // vib defined in Eq. 
(D23). We can now derive Eq. (113) for the transition matrix 
element. This matrix element is within the BO approxima- 
tion given by 



r /0 = ^f(R 1 -..R w )^ (R 1 ...R„), 



(D28) 



N 3 



CMN~ 



1 



2 M k cmn' 

nuc 



f«=2 2 rT r n*(/ r +L AV )a*(^+^) 

k=l irs lM k 



* 1 



t^X 2M k i' r i=\ 
N 



2 ^ k ,{J r +L Ntr )n li p k l +n.c, 



f =y — 



2M 



2 - 



nuc k,n 



P -p" 



(D22) 



where i labels an excited state of the BO Hamiltonian. We 
note that Ar, =r, for i¥=Q and zero otherwise so that we 
only need to consider r, . Through Eq. (D6) we can express 
r /0 in coordinates in the body-fixed frame. Therefore 

| dVr m Qj= j dV<t>* s <t> Qj. (D29) 

The volume element dV has now to be expressed in body- 
fixed frame coordinates. This yields, considering a Jacobian 
consistent with the approximations made in the derivation of 
Eq. (D25), the result 32 
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dV= dTl CMN sm a da dfidydQ, 



(D30) 



where dQ = dQ l - ■ ■ ■ The integral in Eq. (D29) is 

zero if q and / are nonzero. The ground-state center-of-mass 



and rigid-rotor wave functions can then be integrated out to 
yield 

J dVT i0 Qj= ] dQxf(Q)Xo(Q)Qj (D31) 
which proves Eq. (113). 
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